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S ■ Abstract 

• A bifurcation analysis of dune shape transition is made. By use of a reduced model of dune 

morphodynamics, dune skeleton model, we elucidate the transition mechanism between different 
shapes of dunes under unidirectional wind. It was found that the decrease in the total amount of 

■ sand in the system and/or the lateral sand flow shifts the stable state from a straight transverse 

■ dune to wavy transverse dune through a pitchfork bifurcation. A further decrease causes wavy 
transverse dunes to shift into barchans through a Hopf bifurcation. These bifurcation structures 

_ reveal the transition mechanism of dune shapes under unidirectional wind. 
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Sand dunes, considered as the largest granular objects on Earth, have several distinct 
shapeSjfor instance, barchan, transverse, linear, star-shaped, dome-shaped, and parabolic 
dunes [2| . These shapes are governed by two dominant factors; the steadiness of the wind 
direction and the amount of available sand [3]. For example, unidirectional wind generates 
barchans and transverse dunes. The former, crescent-shaped isolated dunes, are formed in 
dune fields with small amounts of available sand, whereas the latter, extending perpendicular 
to the wind direction, are formed in dune fields with larger amounts of available sand than 
barchan-rich fields. 

Recent dune studies have reflected a significant progress in the quantitative analysis of 
dune morphodynamics. For instance, rescaled water tank experiments have successfully 



been conducted to form distinct dune shapes under controlled conditions 

L 

computational models have been used to reproduce various shapes of dunes yy 



6|. Further, 
12 1 . However, 



a theoretical methodology to explain the basic mechanism behind dune shape formation 
beyond a mere numerical reproduction of the shapes is yet to be established except for a 
few studies like simple discussion by Parteli et al. [l^. To address this lack of analytical 
methodology, Niiya et al. proposed the Dune skeleton model (hereafter, DS model) consisting 
of coupled ordinary differential equations, each of which represents the dynamics of two- 
dimensional cross sections (hereafter, 2D-CSs) of a three-dimensional dune Q, [l5|. So far, 
the model has successfully reproduced 3 typical shapes of dunes-straight transverse dune, 
wavy transverse dune, and barchan-depending on the amount of available sand and wind 
strength. In this study, using bifurcation analysis, we investigate the shape transition of 
dunes with the change in relevant parameters. 

The DS model covers the formation processes of barchans and transverse dunes, both of 
which are generated under unidirectional steady wind. This model is roughly based on two 
considerations. First, a lateral array of 2D-CSs is set perpendicular to the wind direction 
with constant intervals. Second, a combination of two forms of sand movement-the intra- 
sand movement within each 2D-CS and the inter-sand movement between neighboring 2D- 
CSs-is considered in governing the macroscopic morphodynamics of dunes. 

Considering the observation that 2D-CSs of barchans and transverse dunes very roughly 
show a scale invariant triangular shape, we assume that the angles of their upwind and 
downwind slopes {6 and (f, respectively) are maintained constant through their time evolu- 
tion, irrespective of their size (Fig. H^a)). We also assume that dunes are located on a flat 



FIG. 1. (a) Intra-2D-CS sand flow, (b) Inter-2D-CS sand flow. 



and hard ground. From these assumptions, the horizontal (i.e., wind directional) position 
and the height of each 2D-CS is uniquely determined if only the coordinate (x, h) of its 
crest is given. Moreover, empirical geometrical constants A, B, and C are set to 1/10, 4/5, 
and 1/5, respectively, reflecting the typical 2D-CS profiles of real barchans and transversed 
dunes. 

As mentioned above, sand flow is classified into two forms: (a) the intra-2D-CS flow 
and (b) the inter-2D-CS flow. The intra-2D-CS flow along the upwind slope is uniquely 
determined if the overcrest sand flux q and the incoming sand flux from the windward 
ground /™ are given (Fig. [It^a)). The overcrest sand flux q deposits along the downwind 
slope or escapes to the leeward inter-dune ground. The deposition ratio Tg in the downwind 
slope, termed sand trapping efficiency [2;|, is basically an increase function of height with 
upper limit Te = I. Considering a real-data-based estimation of Te by Momiji et al. jl6l |. 
we roughly assumed the specific form of Te as 



here h is the non-dimensionalized height of dune with the physical unit of meter [l7j. The 
quantity q reflects the over-dune wind strength; here, we assume g to be a constant, inde- 
pendent of the 2D-CS height. In addition, all of the incoming sand flux from the windward 
ground, deposits on the upwind slope. 

The inter-2D-CS flow Ju(i^j)l Jd(j^i) occurs only between the upwind/downwind slopes of 
neighboring 2D-CSs, i and j. The flux is roughly considered as lateral diffusion depending 
on the height difference, though the consideration of the overlap length of slopes causes 
nonlinearity in this inter-2D-CS flux. The specific forms of Ju{i^j) and Jd(j^i) are described 
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where Aw is the lateral interval between neighboring 2D-CSs and is set as Aw = 1 hereafter 



18| . These quantities correspond to the grey areas in Fig. [IJb) multiplied by the diffusion 



coefficients. Here, the upwind and downwind diffusion coefficients {D^ and D^, respectively) 
control the amount of inter-2D-CS sand flow on respective sides of the slopes and reflect the 
over-dune wind strength. 

With consideration of the above intra- and inter-sand flows, the dynamics of the coordi- 
nate (xj, hi) of 2D-CS's crest is given as a system of coupled ordinary differential equations: 
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We also introduce the annihilation rule of the 2D-CSs at the lateral edges of a dune; This 
rule is required to simulate the shrinking process of dunes. This rule is applicable in the 
cases where hi decreases to /ij = or where the overlap of either the upwind or downwind 
slope between the 2D-CS at the edge and its nearest neighbor vanishes. 

In the previous study, numerical simulations of eqs. ( 12a|) and pbj) have been conducted 
with = 1000 2D-CSs. The lateral boundary condition was set to be periodic, whereas the 
wind directional boundary condition was set such that the total amount of escaped sand 
is uniformly redistributed on the windward ground. These simulations showed that three 
typical shapes of dunes-straight transverse dunes, wavy transverse dunes, and barchans- 




FIG. 2. Analysis setup: (a) Two 2D-CSs selected from the transverse dune obtained by simulation: 
first and second 2D-CSs are represented by the solid and dashed lines, respectively, (b) System of 
two 2D-CSs with lateral periodicity. 



are formed depending on the combination of the initial average height of 2D-CSs (in other 
words, the amount of sand), the overcrest sand q, and the upwind diffusion coefficient D^- 

In this study, we elucidate the transition mechanism between different types of dunes. To 
do so, we conduct a bifurcation analysis of the reduced DS model. The reduced DS model 
consists of two 2D-CSs, each of which is assumed to represent a large number of 2D-CSs used 
in the previous simulation (Fig. [2]). The reduced DS model is described by four variables 
(xi, hi,X2, h2), which are the coordinates of the two 2D-CS's crests. 

The lateral boundary condition is set to be periodic, while the wind directional boundary 
condition is set such that sand escaped from the leeward boundary is redistributed uniformly 
at the windward boundary. Because of these boundary conditions, the total amount of sand 
in the system is preserved. Therefore, the total area of the 2D-CSs, 

_ hi + hi _ Hi 

is maintained constant unless the annihilation of a 2D-CS occurs, and we do not treat the 
annihilation process in this study. 

The relation eq. ([3]) is rewritten as /i2(^i) = A/2ifQ — h\. In addition, by defining the 
relative distance y = X2 — Xi between two 2D-CSs, the reduced DS model is further reduced 
into a two- variable system. Moreover, according to the rule for the redistribution of escaping 
sand, the incoming sand fluxes from windward ground fl"^ are expressed as 

fr = \Y.<ii:^-TE{h,)), 1 = 1,2 

5 



where TE{h) is the deposition ratio in the downwind slope given above. In addition, sym- 
metric relations of the inter-2D-CS flow, Ju{i^2) = —Ju{2~>i) and Jd(2-s>i) = —Jd{i^2)i are 
obtained from the definitions of the four variables constituting the two 2D-CSs. Then, the 
dynamics of the reduced DS model are described as follows: 

dy f TE{h2{hi)) _ TE{hi) \ ^ _ ^^^^ ( 1 1 



dt \ h2{hi) hi J \h2{hi) hi 

-2(BJ,(,^,) + CJ„(,^,))(^ + i-), (4a) 

^ = A [q{TE{hi) - 1) + 2(J,(2^i) - J„(i^2)) + fr] ■ (4b) 



In this simplest form of the reduced DS model, we vary 4 parameters, Hq, q, Du, and Dd. 
Hq controls the total amount of sand in the system; q is the overcrest sand, which reflects 
the wind strength in a dune field; and Du and Dd are the upwind and downwind diffusion 
coefficients, respectively, as introduced in eqs. f|Ta|) and f llb[) . D^, similar to q, is influenced 
by the wind strength, although the correlation between Du and q is not discussed here. 

Now, we conduct the analysis of eqs. fHa|) and fl4bp . We obtain two types of fixed points: 
a) fixed points with hi = /i2 and y = 0, referred to as "straight fixed points" and denoted 
as Sp, and b) fixed points with hi ^ h2 and y 0, referred to as "wavy fixed point" and 
denoted as Wp- To avoid the cases requiring the annihilation rule of 2D-CS, we consider the 
inter-2D-CS fiux function (eqs. (fTa|) and (jlbp ) and search for fixed points of eqs. ( Ha|) and 
f l4bp under the limited conditions where the interaction between two 2D-CS is kept non-zero; 
these conditions are given as 

0<hi< V2Ho, -jV2Ho <y< jV2Ho. (5) 

Sp always exists within these ranges; that is, hi = h2 = Hq and y = serve as a fixed 
point of eqs. (Hal) and fl4bp . independent of control parameters q,Du, and D^. However, 
these control parameters largely affect the stability of Sp and the existence of Wp. Here, 
in order to relate the fixed points Sp and Wp to dune shapes obtained by our previous 
simulation of the full DS model with = 1000 2D-CSs as well as to those in the desert field, 
we check the existence of Wp and the stability of 5*^ and Wp. 

I): If Sp is stable, we consider this state as corresponding to the stable straight transverse 
dunes. The phase of the system in this state is referred to as the "ST phase." 
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FIG. 3. Nullclines, solution trajectories, and vector fields of eq. (|4a|) and ()4bp for (a) Hq = 25, 
(b) Hq = 35, (c) Du = 0.06, and (d) = 0.10. The thin and thick solid lines represent y- 
nullcline and /ii-nullcline, respectively. Both short and longer dashed lines represent the solution 
trajectories starting from the different initial values. The vectors are normalized in magnitude 
while maintaining the direction. The unvarying parameters are set as g = 0.5 and = 0.1 in all 
figures. Other parameters are set as ((a), (b)) = 0.1 and ((c), (d)) Hq = 40. 



II) : If Sf becomes unstable and, in turn, Wps emerge as stable fixed points, we consider this 

transition as corresponding to that from the straight transverse dunes to the stable 
wavy transverse dunes. The phase of the system in this state is referred to as the "WT 
phase." 

III) : If Wp becomes unstable with Sp already unstable, we consider this transition as 

corresponding to that from wavy transverse dunes to barchans. The phase of the 
system in this state is referred to as the "B phase." 

In the B phase, no stable fixed point is obtained. Correspondingly, a single barchan is not 
stable in the DS model simulation and gradually shrinks with time. Note that the later- 
ally periodic boundary condition of the present model should be reconsidered to correctly 
describe the dynamics of barchans. 

Now, we study the phase change transition between ST and WT phases and that between 
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FIG. 4. Bifurcation diagram of fixed points obtained by varying one control parameter: (a) Hq 
and (b) Du- The solid and dashed lines denote stable and unstable fixed points, respectively. B, W, 
and S denote the ranges of B phase, WT phase, and ST phase, corresponding to the dune shapes, 
respectively. The unvarying parameters {q, D^) that are common between (a) and (b) are the same 
as the values in Fig. [3l whereas other parameters are set as (a) Du = 0.1 and (b) Hq = 40. 



WT and B phases using the bifurcation analysis. First, we consider the nullchnes of eqs. 
f Ha|) and fl4b|) as a two-dimensional dynamical system, where a pair of control parameters 
are fixed as g = 0.5 and = 0.1 and the remaining two parameters, Hq and D^, are varied 
independently. In Fig. |3l the nullclines in the y-hi phase space are drawn at typical values 
of Hq or Du- They suggest that the decrease in Hq and/or Du plays a significant role in the 
emergence of Wp. Here, the results are unchanged even if the two 2D-CSs are exchanged; 
that is, below a critical value of Hq or Du, two points of Wp simultaneously emerge and are 
symmetrical around Sp (Fig. Et^a), (c)). 

Next, in order to determine the bifurcation structure and stability of fixed points, we 
search for fixed points by continuously varying the parameters Hq and Du within the range 
satisfying relation eq. ([5]) and conduct the linear stability analysis (Fig. H]). The bifurcation 
diagrams clearly show the transitions between three phases: ST, WT, and B phases. Figure 
m^a) and m^b) shows that the decrease in Hq and/or Du destabilizes Sp to shift the corre- 
sponding phase of the system from ST to WT and further decrease in the same destabilizes 
Wp to shift the system into B phase. It should be noted that B phase in Fig. Hl^b) occurs 
when Wp goes out of the range of eq. ([5]). 

Moreover, in Fig. ID^a) andHJ^b), we focus on two bifurcation points in the respective 
figures: (i) first bifurcation points of switching from ST phase to WT phase and (ii) second 
bifurcation points of switching from WT phase to B phase. The type of bifurcation at the 
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first bifurcation points is tlie pitclifork bifurcation because the stable Sp becomes unstable 
and two symmetric stable Wp points arise with the decrease in Hq and Du- Further decrease 
in Hq and Du destabilizes Wp, however, it does not cause the emergence of new fixed points. 
Beyond these second bifurcation points, the eigenvalues of the Jacobian matrix at Wp contain 
imaginary parts from the real phase space. This type of bifurcation at the second bifurcation 
points is the Hopf bifurcation. 

Finally, we classify each of the two sets of parameter spaces [Hq, q) and [D^, q) into three 
phases using the bifurcation analysis discussed above. These phase diagrams qualitatively 
correspond to our previous numerical simulation results of the DS model [l^, [3] • Namely, 
an increase in the amount of sand and inter-2D-CSs flow enhances the stability of the 
transverse dune, whereas the decrease in intra-2D-CS flow destabilizes its shape to enforce 
the deformation to a barchan. Moreover, these shape transition reproduce the change in 
real dune shapes depending on the amount of available sand. 

In this study, we dealt with dune morpho dynamics from a dynamical system point of 
view. A simple model called the Dune Skeleton model was further simplified to a two- 
variable ordinary differential equation named the ^^reduced Dune Skeleton model." This 
reduced DS model enabled the bifurcation analysis for the transition between different dune 
shapes observed under a unidirectional wind. The results qualitatively corresponded to the 
previous observations of real dunes, water tank experiments, and simulation models. But 
the mere careful comparison to recent study which are made under similar but not identical 
condition is next issue. 

This work was partially supported by the 21st Century COE Program, "Toward a new 
basic science depth and synthesis." 
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